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The hierarchical equations of motion theory for Drude dissipation is optimized, with a convenient 
convergence criterion proposed in advance of numerical propagations. The theoretical construction 
is on basis of a Pade spectrum decomposition that has been qualified to be the best sum-over-poles 
scheme for quantum distribution function. The resulting hierarchical dynamics under the apriori 
convergence criterion are exemplified with a benchmark spin-boson system, and also the transient 
absorption and two-dimensional spectroscopy of a model exciton dimer system. 
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I. INTRODUCTION 

Quantum dissipation plays a pivotal role in many 
chemical, physical, and biological problems ii"— Recent 
experiments on excitation energy transfer in photosyn- 
thesis systems^i^ show clearly the breakdown of conven- 
tional Markovian and perturbative quantum dissipation 
theoriesi^ The protein environment is nano-structured, 
with the time scale of modulation comparable to that of 
the excitation energy transfer. Also the protein-pigment 
coupling strength is about that between pigments. Ap- 
parently, one needs a nonperturbative, non-Markovian, 
and also numerically implementable quantum dissipation 
theory. The demand is the same on the study of heat- 
ing generation in quantum transport through mcsoscopic 
systems2r2i and protection of qubits with the aid of pho- 
tonic crystal environment J^"— 

This work focuses on the hierarchical equations of 
motion (HEOM) approachfi^"— We request the best 
HEOM construction, exemplified with Drude dissipa- 
tion, together with a convenient criterion to estimate 
its convergency before propagations for general sys- 
tems at any finite temperature. HEOM approach was 
originally proposed in 1989 by Tanimura and Kubo 
for semiclassical dissipation)^ Formally exact HEOM 



formalisi 
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for Gaussian dissipation in general, in- 



cluding its second quantization^ is now well established, 
with the aid of proper environment spectrum decompo- 
sition schemes. As a numerically efficient alternative to 
path integral influence functionalf^i^i^ HEOM has been 
applied to such as electron transfer f2ir— nonlinear optical 
spectroscopy)^"— and transient quantum transport i^"— 

To have an explicit HEOM construction one should ex- 
ploit certain basis set spanning over the stochastic bath 
space. This is equivalent to the choice of certain sum- 
over-pole (SOP) scheme that decomposes individual bath 
correlation function into its multiple-timescale memory 
spectrum components. The conventional scheme is the 
Matsubara expansion of quantum distribution function 
(i.e., the Bose/Fermi function) involved in the bath corre- 
lation functions.— — However, Matsubara expansion is no- 
torious for slow convergence. We have recently proposed 
three Fade spectrum decomposition (FSD) schcmes^i^ 



be the candidates for the best SOF method. Mathemat- 
ically, these three FSDs of Bose/Fermi function exploit 
the [7V-l/iV], [N/N], and [N+l/N] Fade approximants, 
respectively. The resulting HEOM dynamics have been 
demonstrated in context of transient quantum transport 
through a double quantum dots system and population 
transfer in a spin-boson systemi^ 

Three closely related issues arise in the choice of statis- 
tical environment basis set (or SOF scheme) for efficient 
HEOM construction. The first one is to identify the best 
basis set for spanning over the stochastic bath space. The 
second one concerns the minimum basis set. It requests 
to have not only the smallest number of decomposition 
terms, but also a priori accuracy control criterion on the 
resulting HEOM dynamics for any given system under 
bath influence at finite temperature. The third issue is 
about the possibility of at least partial inclusion of the 
off-basis-set residue effect on the HEOM dynamics. 

In this work, we construct an optimized HEOM the- 
ory for Drude dissipation. We identity that [N/N] FSD 
serves as the best and the minimum Drude dissipation ba- 
sis set. It leads naturally to an off-basis-set white noise 
residue (WNR) term, together with a simple accuracy 
control criterion on the resulting HEOM dynamics. The 
present paper is a generalization of Ref. |25| and Ref. |26| . 
where A^ = and A^ = 1 cases were analyzed, respec- 
tively. We summarize the FSD-based HEOM formalism 
with the Drude bath in SecHH followed by proposing the 
accuracy control criterion. Numerical demonstrations are 
carried out in Sec. lIIII Included there are a benchmark 
spin-boson dynamics, studied before by Thoss, Wang, 
and Miller;^ and the nonlinear optical signals of a model 
dimer system. Finally we conclude the paper in Sec. lIVI 



II. FORMALISM 

A. Optimal hierarchy for Drude dissipation 

In this subsection, we exploit the [N/N] FSD scheme^ 
to construct HEOM for Drude dissipation cases. HEOM 
has the following generic fornifSfl 

Pn(t) = -[iCit) + 7„ + 6n,]p„{t) + pV (t) + p^' (t). (1) 



It describes how an n*'^-tier ADO pn depends on its asso- 
ciated (n± l)*^-tier ADOs in p,^**. The ADO's index n is 
in general a coUection of indices; i.e., n = {rii, • • • ,rix}, 
with rik > for bosonic bath. Here ni + ■ ■ ■ + tik = n 
for pn = Pm.-- .UK being called an n"^-tier ADO. The re- 
duced system density operator p(t) = trBptotai(0 — Po(i) 
is just the zeroth-tier ADO. In Eq. ([1]), the reduced sys- 
tem Liouvillian C{t) ■ = [H{t)^ ■ ] can be time dependent, 
e.g., in the presence of driving fields. Throughout of this 
paper, we set h = 1 and /3 = l/ikgT), with kg being the 
Boltzmann constant and T the temperature. 

The specific HEOM construction, including the ADO 
labeling index n = {ni, • • • , nx}, depends on the way of 
decomposing bath correlation function into its memory 
spectrum components. For clarity, let the system-bath 
interaction be H'{t) = —QF^{t), with Q and F-s{t) be- 
ing operators in the reduced system and the stochastic 
bath subspaces, respectively. The influence of bath on 
system is completely determined by the correlation func- 
tion C{t) = {F^{t)F-a({)))^. It is in turn related to the 
bath spectral density J{ljj) via the fluctuation-dissipation 
theorem (FDT)i^-l 
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To have a HEOM construction,— we expand C{t) in a 
finite exponential series, on the basis of certain SOP 
scheme, together with the Cauchy residue theorem of 
contour integration. 

In this work we focus on the Drudc model, 



J(^) 
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It has only one pole, z ~ —17 = —470, hi the lower-half 
plane. The exponential series of bath correlation function 
assumes then 



TV 
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SCN{t). 



(4) 
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In general the ofF-basis-set residue SCpf{t) 7^ 0, as one 
exploits only finite N poles for Bose function. The con- 
ventional scheme is the Matsubara expansion which is 
however notorious for slow convergence. ^^ For Drude dis- 
sipation it has been suggestedS^ that [N/N] PSD be the 
best SOP of Bose function. It is accurate up to 0{x'^^^^) 
in the order of a; = /3uj and reads^ 
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with 



Rn = 



4{N + l){2N + 3) 



'Rnx, (5) 



(6) 



The PSD poles and coefficients, {(,k,Vk',k — !,■■■ ,N}, 
are all positive and can be evaluated with high precision 
via the eigenvalues of real symmetric triangle matrices j^ 



The corresponding exponential series expansion of 
Eq. ^ can now be obtained via the standard contour 
integration technique. We obtain (setting jk = £.k/P) 
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and 



SC'Nit) « 4:Xf3-fRN6{t) = 2AAr(5(<). 



(7) 



(8) 



Note that {cfc; fc = 1, • • • , N} are all real. On the other 
hand, Cd associating with the Drude exponent 7d = 7 is 
complex and is evaluated by using the [N/N] Bose func- 
tion /'^/^'(x), as given by the last expression of Eq. ([5]). 
The SOP scheme resembles a statistical environment 
basis set for the HEOM construction. It dictates not only 
the exponential expansion of C{t) as Eq. ^, but also 



HEOM. The ADO reads now p„ = /9„d,„i, 



The 



only approximation involved is the WNR treatment of 
the off-basis-set bCnii:) by Eq. ^. It results in the WNR 
of (57?.n in Eq. (H)) without further approximation! ^^'^° 



^TlnPxx = '^n[Q, [Q,/On] 



(9) 



The damping parameter in Eq. ([T]) collects all relevant 
exponents 
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with p ± being the associated {n ± l)*^-tier ADO, re- 

k 

spectively. The labeling index n^ differs from n only by 
changing the specified Uk to rife ± 1. All ADOs here are 
dimensionless and scaled properly for the efficient HEOM 
propagator with the recently developed on-the-fly filter- 
ing algorithm4i Numerically it also automatically trun- 
cates the level of hierarchy4i 



B. Accuracy control criterions 

As mentioned earlier the only approximation involved 
is the WNR treatment of the off-basis-set SCN{t) by 
Eq. ([8|) . Its validity dictates therefore the accuracy of 
the resulted HEOM. The exact residue 5Cisi{t), which 
is a real and even function, is defined via Eq. ([3]), to- 
gether with Eq. ([7]) for the Drude model. Its spectrum 



6Cn{i^) = 5 /<iie*'^*(5CAr(i) is a symmetric and bell- 
shaped function, being positive and monotonically de- 
creasing in w > from SCn{^ = 0) = Ajv to 5Cn{oj — J- 
oo) = 0, where Aat = 2Xj3jRn was defined in Eq. ([5]). 
The half- width-at-half- maximum T^('y) that character- 
izes the inverse time scale of 5CM{t) is determined via 
5Cn{^)\uj=Tn(i) ~ ^Af/2- It is found that r7v(7) can be 
well approximated by (within 0.5% of relative error for 
iV < 16 as tested) 



^n{i) ~ ^ 



VN 



(/37)2 + 0.34r2 
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(12) 



where tn = 1/{2Rn) = 2(iV + \){2N + 3) [cf. Eq. ©]. 

Figure [T] depicts the residue spectrum 5Cn{uj), plot- 
ted in terms of SCn{(^)/^n versus uj/Tis[('j) for some 
selected values of {N, P'^}. Used here for the x-axis scal- 
ing is the approximated Fjv (7) value of Eq. P^ . Thus 
this figure shows also the excellent quality of Eq. P^ . 

To validate the 5- function approximation of 5CM{t) 
[Eq. (|S])], we examine the Kubo's motional narrowing line 
shape parameter )^2ii^ 



Kiv(7, A) = ^Tn{i)I^n = yJrNTN{l)l{P\l) ■ (13) 

The evaluated I3Tn{'^) and Sat (7) = (/3A)^/^K7v(7, A), as 
functions of /37, are depicted in Fig.[51Ja) and (b), respec- 
tively. 

The accuracy control criterion on the HEOM in 
Sec. Ill A] comprises therefore the conditions under which 
5CN{t) and its effect on the reduced system dynamics 
can be treated as Markovian white noise. Apparently, 
it is valid when FAr(7) 3> £7^ and KAr(7, A) 3> 1, with 
r2s denoting the characteristic frequency of system. We 
have foun d^^'^^ that the HEOM dynamics assumes nu- 
merically accurate (if not exact) when 



i{Fw(7)M,Kw(7,A)}>5, 



(14) 



while semi-quantitative when 

2 < min{Fjv(7)M, ^^^(7, A)} < 5. (15) 

The above accuracy control criterions facilitate the choice 
of minimum N for the desired quality of HEOM dynam- 
ics; see Refs. [2^ and [2a for the special cases of A^ = 
and 1, respectively. 

III. NUMERICAL DEMONSTRATIONS 

A. Spin-boson dynamics 

To demonstrate the efficiency of HEOM and also the 
proposed accuracy control criterions, consider a bench- 
mark spin-boson system, studied before by Thoss, Wang, 
and Miller via their numerically exact self-consistent hy- 
brid approach)^ The system Hamiltonian is H = eaz + 
Vax, with the dissipative mode of Q — az, subject to 



Drude dissipation. Here ctx and ctz are Pauli matrices. 
The Rabi frequency of the system is fig = 2Ve^ + V^. 
We choose the most challenging parameters used in Ref. 
EO, i.e., the figure 8 there, with e/V = 1, X/V = 0.25, 
7/y = 5, and ^V = 50. 

Figure [3] depicts the reduced system density matrix 
evolution, as evaluated from the HEOM at each spec- 
ified environment basis set size N. It converges when 
N > 10 (by eye inspection). The accuracy control pa- 
rameters {TN/rts,KN} are {2.6,3.6}Ar=4, {4.7,8.5}Ar=8; 
{6.3, 12.0}jv=io, and {8.4, 16.3}Ar=i2, respectively. The 
CPU time listed in Fig.|3] is based on a single Intel(R) 
Xeon(R) processor@3.00GHz. The HEOM propagation 
uses the fourth-order Runge-Kutta method with time- 
step of 0.001/y, together with the on-the-fly filtering 
algorithm^ with error tolerance of 5 x 10~^. Evidently 
the numerically accurate criterion (J14p and the semi- 
quantitative one (J15p are both verified. Moreover, the 
[N/N] PSD-based HEOM is found to be the best among 
all schemes we tested. 



B. Nonlinear spectroscopy of a model exciton 
dimer system 

We now examine the accuracy control over the [N/N] 
PSD-based HEOM in evaluating nonlinear spectroscopic 
signals. We will verify again the proposed accuracy con- 
trol criterions in Eqs. p^ and ((TSj) . respectively. 

We first re-evaluate the figures 6 and 7 of Ref. [2y, due 
to a careless mistake in the excitation field configura- 
tion used therei^ The correct signals (in rotating wave 
approximation) are now depicted in Figs.|4] and [5l re- 
spectively. Demonstrated are different components of 
the dispersed transient absorption coefficient a(w,td) of 
a model exciton dimer system at two representing tem- 
peratures. Here, tj, denotes the probe delay time with re- 
spect to the pump pulse. Denote also Aw = uj — e, where 
e = ei = £2 is the on-site excitonic energy. The dimer sys- 
tem parameters^ are V = — 400cm~^ for exciton trans- 
fer, U = 200 cm""'^ for exciton Coulomb interaction, and 
M2z/miz — 0-35 for the dimer transition dipoles orienta- 
tion asymmetry. The characteristic frequency of system 
is fis = \/(ei - £2)2 -h 4^2 =: 2V. The on-site Drude 
fluctuation parameters are 7 = 71 = 72 = 600 cm~^ and 
A = Ai = A2 = 600cm-i. 

The pump field is a transform-limited Gaussian pulse 
with 50 fs at the full width at half maximum, centered 
at w = e; i.e.. Aw = w — e = 0, as indicated by the 
arrow in panel (a) of each Fig.S] and Fig.[S] The peak 
intensity is of /ii^E'max = 100 cm~^. At td = lOOfs the 
pump-transferred occupations in the single on-site exci- 
ton |1) and |2) states, and the doublc-exciton |/) state 
are 6.0%, 5.9%, and 1.7%, respectively, at T = 298 K; 
while they are 6.3%, 6.1%, and 1.9%, respectively, at 
T = 77K. The probe field assumes in the weak and im- 
pulsive limit. Apparently, the dips appearing in the non- 
linear absorptive a^'^(w, td) components in the (b) panels 



arise mainly from the excite-statc absorption to the dou- 
bly excited state. The linear emission signal [dash-curve 
in (a)] involves the transition from the lower lying sin- 
gle exciton eigenstate to ground state, without involving 
the the double-exciton |/) state. However, as the mod- 
erately strong pump field is used, the nonlinear emissive 
a^^(uj,td) component [(c) panels] contains (in the blue 
side) also the contribution of the [/) state emission. 

Figure [6] shows the two-dimensional (2D) spectroscopy 
S{u!3, td,iOi) of the same dimer system at T = 298 K. As 
inferred from Fig.H the [1/1] PSD-based HEOM is suffi- 
cient to give the numerically accurate results here. Both 
pump and probe fields are now operated in the weak and 
impulsive limit. As a result, S{lo3, td, uji) resolves simply 
both the excitation and detection frequencies, wi and 
L03, of the third-order optical response function, where i2 
is just the delay time td of detection.-^ For demonstra- 
tion we examine again the absorptive (upper panels) and 
the emissive (middle panels) components of the experi- 
mentally measurable S{uj3,td,uJi) (bottom panels) that 
amounts to the kj + kjj signal4i Evidently the dips ap- 
pearing in the absorptive components are due to the ex- 
cited state absorption, while the 2D emissive pathways 
do not include the |/)-state emission, as the pump field 
is now in the weak response regime. 

Figure [7] depicts the Awi = slices of the absorptive 
and emissive panels of Fig. [5] It thus resembles the im- 
pulsive pump counterpart of Fig.|4l except for the emis- 
sive contribution from the doubly excited state. Un- 
like the impulsive limit that considers only sequential 



processes;^ the evaluation with finite pulse duration in- 
volves also coherent processes that cannot be neglected 
at the short time region {td ^ here) 4^"— 



IV. SUMMARY 

In summary, we have constructed the [N/N] PSD- 
based HEOM, which is qualified to be the best hierar- 
chical theory for Drude dissipation. The present work 
generalizes the previous [0/0] and [1/1] PSD-based hier- 
archical constructions i2^^ The proposed accuracy con- 
trol critcrions [Eqs. (|14p and (|15p ] are confirmed via not 
only the reduced system density matrix dynamics but 
also nonlinear optical spectroscopy calculations. No ex- 
pensive convcrgcncy check would thus be needed for the 
HEOM dynamics of complex open systems. HEOM for 
non-Drude environments that should be optimized case 
by case together with accuracy control critcrions will be 
considered elsewhere. 
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FIG. 2: /3riv(7) and KN{-f) = (/3A)^/Viv(7, A) versus /Sj. 
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FIG. 1: The deviation spectrum function SCn{u}) of Drude 
bath, plotted in terms of JCjv(i^)/Ajv versus a;/rAr(7) for 
some selected values of {A'^, P^}, where Am is given via Eq. ((8| 
and rjv(7) is by the approximant of Eq. (|12p. 
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FIG. 3: Evolution of the reduced spin system density ma- 
trix elements, with the same parameters as Fig. 8 of Ref. mj : 
i.e., e/V = 1, X/V = 0.25, j/V = 5, and jSV = 50. The 
accuracy control parameters are {rjv/Jla, kjv} — {2.6,3.6} 
for iV = 4 and {4.7,8.5} for A^ = 8, as indicated. The 
converged dynamics are identical to that of iV = 10 whose 
{rjvM,Kiv} = {6.3,12.0}. 






3 

■J 
z < 



-2000 




2000 



FIG. 4: Dispersed transient absorption coefficient signals 
of the model dimer system (see text): (a) Linear absorp- 
tion (A) and emission (E) signals; (b) Nonlinear absorp- 
tive ci^(ij^^td) component; (c) Nonlinear emissive (J^iy^^td) 
component. The pump field is a transform-limited 50 fs- 
pulse of finite intensity (see text), centered at oj = e; i.e., 
Atj = (X) — e = 0, as indicated by the arrow in panel (a). At 
T = 298 K, A^ = 1 is in the numerically accurate range, with 
{rjv/f2s, Kiv} = {8.3, 8.7} as indicated. 
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FIG. 5: Same as Fig.H but at T = 77 K, where A = 1 is in 
the semiquantitative range, with {r]v/ns,Kjv} = {2.4,2.4}, 
as indicated. 
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FIG. 6: Two-dimensional spectroscopy S{iJs,td,ijJi) of the 
same dimer system of Fig.|3] at T = 298 K. Both pump and 
probe fields are weak and in the impulsive limit. 
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FIG. 7: The Aij;i — slices of the absorptive and emissive 
panels of Fig.[6l This figure resembles the impulsive pump 
counterpart of Fig.|4l 



